rm(list = ls())
source("Accepted Paper Replication Archive/code/utils.R")


data <- read_csv(file =  "Accepted Paper Replication Archive/clean_data/pilot_analysis_data.csv")
head(data)
data <- data[-c(1:2), ] ## remove rows of labels
dim(data)

table(duplicated(data$prolific_ID))

names(data)

ID1 <- unique(data$ResponseId)


###
### DEMOS
###

table(data$age, exclude = NULL)
data$age2 <- as.numeric(as.character(data$age))
data$age2[data$age == "100 or over"] <- 100
summary(data$age2)
dim(data)
data <- data[data$age != "Under 18", ] ## drop respondents under 18
dim(data)

table(data$gender, exclude = NULL)
data$female <- NA
data$female[data$gender == "Female"] <- 1
data$female[data$gender == "Male"] <- 0
data$female[data$gender == "Other"] <- NA ## too small a group to analyze
table(data$female, exclude = NULL)

## examine 'other' responses to gender
table(data$gender_4_TEXT[data$gender == "Other"])

## dichotomous party measure
table(data$party, exclude = NULL)
data$party2 <- NA
data$party2[data$party == "Democrat"] <- "dem"
data$party2[data$party == "Republican"] <- "rep"
data$party2[data$closer == "Democratic Party"] <- "dem"
data$party2[data$closer == "Republican Party"] <- "rep"
data$party2[data$closer == "Neither"] <- "ind"
table(data$party2, exclude = NULL)


###
## TREATMENTS
###
# Dose response experiment
table(data$exp_1_gender) ## gender % emphasized
data$e1_gender <- NA
data$e1_gender[data$exp_1_gender == "men"] <- "men"
data$e1_gender[data$exp_1_gender == "women"] <- "women"
data$e1_gender <- factor(data$e1_gender, levels = c("men", "women"))
class(data$e1_gender)
table(data$e1_gender, exclude = NULL)

table(data$pct_dose, exclude = NULL)
data$pct_dose <- as.numeric(as.character(data$pct_dose))
data$pct_dose <- factor(data$pct_dose, levels = c("50", "20", "30", "40", "60", "70", "80"))
table(data$pct_dose)

table(data$pct_dose)
data$pct_dose2 <- NA
data$pct_dose2[(data$pct_dose == "20" & data$e1_gender == "men") == T |
  (data$pct_dose == "80" & data$e1_gender == "women") == T] <- "80_20"
data$pct_dose2[(data$pct_dose == "80" & data$e1_gender == "men") == T |
  (data$pct_dose == "20" & data$e1_gender == "women") == T] <- "20_80"

data$pct_dose2[(data$pct_dose == "30" & data$e1_gender == "men") == T |
  (data$pct_dose == "70" & data$e1_gender == "women") == T] <- "70_30"
data$pct_dose2[(data$pct_dose == "70" & data$e1_gender == "men") == T |
  (data$pct_dose == "30" & data$e1_gender == "women") == T] <- "30_70"

data$pct_dose2[(data$pct_dose == "40" & data$e1_gender == "men") == T |
  (data$pct_dose == "60" & data$e1_gender == "women") == T] <- "60_40"
data$pct_dose2[(data$pct_dose == "60" & data$e1_gender == "men") == T |
  (data$pct_dose == "40" & data$e1_gender == "women") == T] <- "40_60"

data$pct_dose2[(data$pct_dose == "50" & data$e1_gender == "men") == T |
  (data$pct_dose == "50" & data$e1_gender == "women") == T] <- "50_50"
table(data$pct_dose2)
data$pct_dose2 <- factor(data$pct_dose2, levels = c(
  "20_80", "30_70", "40_60", "50_50",
  "60_40", "70_30", "80_20"
))
table(data$rank)
data$rank <- factor(data$rank, levels = c("jobs", "top jobs"))
table(data$rank)

table(data$exp1_agency) ## agency
data$e1_agency <- NA
data$e1_agency[data$exp1_agency == "Department of Health and Human Services"] <- "hhs"
data$e1_agency[data$exp1_agency == "Department of Defense"] <- "defense"
data$e1_agency[data$exp1_agency == "Department of Education"] <- "education"
data$e1_agency[data$exp1_agency == "Department of Treasury"] <- "treasury"
table(data$e1_agency, exclude = NULL)


## Exp 2
## Nominee experiment
table(data$exp2_agency) ## agency treatment
data$e2_agency <- NA
data$e2_agency[data$exp2_agency == "Department of Health and Human Services"] <- "hhs"
data$e2_agency[data$exp2_agency == "Department of Defense"] <- "defense"
data$e2_agency[data$exp2_agency == "Department of Education"] <- "education"
data$e2_agency[data$exp2_agency == "Department of Treasury"] <- "treasury"
table(data$e2_agency, exclude = NULL)

table(data$exp2_gender) ## gender treatment
data$e2_gender <- NA
data$e2_gender[data$exp2_gender == "man"] <- "male"
data$e2_gender[data$exp2_gender == "woman"] <- "female"
table(data$e2_gender, exclude = NULL)

### Exp 3
### Chart experiment
table(data$FL_174_DO) ### chart displayed
data$e3_chart <- NA
data$e3_chart[data$FL_174_DO == "chart1"] <- "top_jobs_levels"
data$e3_chart[data$FL_174_DO == "chart2"] <- "top_jobs_change_zoom"
data$e3_chart[data$FL_174_DO == "chart2b"] <- "top_jobs_change"
data$e3_chart[data$FL_174_DO == "chart3"] <- "top_rf_levels"
data$e3_chart[data$FL_174_DO == "chart4"] <- "top_rf_US_levels"
data$e3_chart <- factor(data$e3_chart, levels = c(
  "top_jobs_change", "top_jobs_change_zoom",
  "top_jobs_levels", "top_rf_levels",
  "top_rf_US_levels"
))
table(data$e3_chart, exclude = NULL)



### OUTCOMES

# Dose response experiment
## confidence in agency's ability to fulfill mission
table(data$ability_dept)
data$e1_agency_confidence <- NA
data$e1_agency_confidence[data$ability_dept == "No confidence at all"] <- 1
data$e1_agency_confidence[data$ability_dept == "Not too much confidence"] <- 2
data$e1_agency_confidence[data$ability_dept == "Some confidence"] <- 3
data$e1_agency_confidence[data$ability_dept == "A lot of confidence"] <- 4
table(data$e1_agency_confidence, exclude = NULL)
data$e1_agency_confidence <- recode_0_1(data$e1_agency_confidence)
table(data$e1_agency_confidence, exclude = NULL)


## confidence agency will rep interests of women
table(data$rep_women)
data$e1_rep_women <- NA
data$e1_rep_women[data$rep_women == "No confidence at all"] <- 1
data$e1_rep_women[data$rep_women == "Not too much confidence"] <- 2
data$e1_rep_women[data$rep_women == "Some confidence"] <- 3
data$e1_rep_women[data$rep_women == "A lot of confidence"] <- 4
table(data$e1_rep_women, exclude = NULL)
data$e1_rep_women <- recode_0_1(data$e1_rep_women)
table(data$e1_rep_women, exclude = NULL)


# Nominee experiment
## ability of nominee to effectively lead agency
table(data$exp2_dv1)
data$e2_nominee_confidence <- NA
data$e2_nominee_confidence[data$exp2_dv1 == "No confidence at all"] <- 1
data$e2_nominee_confidence[data$exp2_dv1 == "Not too much confidence"] <- 2
data$e2_nominee_confidence[data$exp2_dv1 == "Some confidence"] <- 3
data$e2_nominee_confidence[data$exp2_dv1 == "A lot of confidence"] <- 4
table(data$e2_nominee_confidence, exclude = NULL)
data$e2_nominee_confidence <- recode_0_1(data$e2_nominee_confidence)
table(data$e2_nominee_confidence, exclude = NULL)


## president's ability to effectively staff government
table(data$exp2_dv2)
data$e2_pres_confidence <- NA
data$e2_pres_confidence[data$exp2_dv2 == "No confidence at all"] <- 1
data$e2_pres_confidence[data$exp2_dv2 == "Not too much confidence"] <- 2
data$e2_pres_confidence[data$exp2_dv2 == "Some confidence"] <- 3
data$e2_pres_confidence[data$exp2_dv2 == "A lot of confidence"] <- 4
table(data$e2_pres_confidence, exclude = NULL)
data$e2_pres_confidence <- recode_0_1(data$e2_pres_confidence)
table(data$e2_pres_confidence, exclude = NULL)


## women have adequate rep in government
table(data$exp3_dv1)
table(data$exp3_dv11)
table(data$exp3_dv1b)
table(data$exp3_dv1c)
table(data$exp3_dv1d)

data$e3_women_rep <- NA
data$e3_women_rep[data$exp3_dv1 == "Much too little representation"] <- 1
data$e3_women_rep[data$exp3_dv11 == "Much too little representation"] <- 1
data$e3_women_rep[data$exp3_dv1b == "Much too little representation"] <- 1
data$e3_women_rep[data$exp3_dv1c == "Much too little representation"] <- 1
data$e3_women_rep[data$exp3_dv1d == "Much too little representation"] <- 1

data$e3_women_rep[data$exp3_dv1 == "Too little representation"] <- 2
data$e3_women_rep[data$exp3_dv11 == "Too little representation"] <- 2
data$e3_women_rep[data$exp3_dv1b == "Too little representation"] <- 2
data$e3_women_rep[data$exp3_dv1c == "Too little representation"] <- 2
data$e3_women_rep[data$exp3_dv1d == "Too little representation"] <- 2

data$e3_women_rep[data$exp3_dv1 == "About the right amount of representation"] <- 3
data$e3_women_rep[data$exp3_dv11 == "About the right amount of representation"] <- 3
data$e3_women_rep[data$exp3_dv1b == "About the right amount of representation"] <- 3
data$e3_women_rep[data$exp3_dv1c == "About the right amount of representation"] <- 3
data$e3_women_rep[data$exp3_dv1d == "About the right amount of representation"] <- 3

data$e3_women_rep[data$exp3_dv1 == "Too much representation"] <- 4
data$e3_women_rep[data$exp3_dv11 == "Too much representation"] <- 4
data$e3_women_rep[data$exp3_dv1b == "Too much representation"] <- 4
data$e3_women_rep[data$exp3_dv1c == "Too much representation"] <- 4
data$e3_women_rep[data$exp3_dv1d == "Too much representation"] <- 4

data$e3_women_rep[data$exp3_dv1 == "Much too much representation"] <- 5
data$e3_women_rep[data$exp3_dv11 == "Much too much representation"] <- 5
data$e3_women_rep[data$exp3_dv1b == "Much too much representation"] <- 5
data$e3_women_rep[data$exp3_dv1c == "Much too much representation"] <- 5
data$e3_women_rep[data$exp3_dv1d == "Much too much representation"] <- 5
table(data$e3_women_rep, exclude = NULL)
data$e3_women_rep <- recode_0_1(data$e3_women_rep)
table(data$e3_women_rep, exclude = NULL)


## US Gov

table(data$exp3_dv3[1:162], exclude = NULL)
table(data$exp3_dv32)
table(data$exp3_dv3b)
table(data$exp3_dv3c)
table(data$exp3_dv3d)

data$e3_serve_US <- NA
data$e3_serve_US[data$exp3_dv3 == ""] <- NA
data$e3_serve_US[data$exp3_dv32 == ""] <- NA
data$e3_serve_US[data$exp3_dv3b == ""] <- NA
data$e3_serve_US[data$exp3_dv3c == ""] <- NA
data$e3_serve_US[data$exp3_dv3d == ""] <- NA
table(data$e3_serve_US, exclude = NULL)


data$e3_serve_US[data$exp3_dv3 == "No confidence at all"] <- 1
data$e3_serve_US[data$exp3_dv32 == "No confidence at all"] <- 1
data$e3_serve_US[data$exp3_dv3b == "No confidence at all"] <- 1
data$e3_serve_US[data$exp3_dv3c == "No confidence at all"] <- 1
data$e3_serve_US[data$exp3_dv3d == "No confidence at all"] <- 1
table(data$e3_serve_US, exclude = NULL)


data$e3_serve_US[data$exp3_dv3 == "Not too much confidence"] <- 2
data$e3_serve_US[data$exp3_dv32 == "Not too much confidence"] <- 2
data$e3_serve_US[data$exp3_dv3b == "Not too much confidence"] <- 2
data$e3_serve_US[data$exp3_dv3c == "Not too much confidence"] <- 2
data$e3_serve_US[data$exp3_dv3d == "Not too much confidence"] <- 2
table(data$e3_serve_US, exclude = NULL)


data$e3_serve_US[data$exp3_dv3 == "Some confidence"] <- 3
data$e3_serve_US[data$exp3_dv32 == "Some confidence"] <- 3
data$e3_serve_US[data$exp3_dv3b == "Some confidence"] <- 3
data$e3_serve_US[data$exp3_dv3c == "Some confidence"] <- 3
data$e3_serve_US[data$exp3_dv3d == "Some confidence"] <- 3
table(data$e3_serve_US, exclude = NULL)


data$e3_serve_US[data$exp3_dv3 == "A lot of confidence"] <- 4
data$e3_serve_US[data$exp3_dv32 == "A lot of confidence"] <- 4
data$e3_serve_US[data$exp3_dv3b == "A lot of confidence"] <- 4
data$e3_serve_US[data$exp3_dv3c == "A lot of confidence"] <- 4
data$e3_serve_US[data$exp3_dv3d == "A lot of confidence"] <- 4
table(data$e3_serve_US, exclude = NULL)


table(data$e3_serve_US, exclude = NULL)
data$e3_serve_US <- recode_0_1(data$e3_serve_US)
table(data$e3_serve_US, exclude = NULL)

### order of experiments

## FL_133 = dose response
## FL_178 = job candidate
## FL_169 = charts
data$order_exp <- NA
data$order_exp[data$FL_134_DO == "FL_133|FL_169|FL_178" | data$FL_134_DO == "FL_133|FL_178|FL_169"] <- "dose"
data$order_exp[data$FL_134_DO == "FL_178|FL_133|FL_169" | data$FL_134_DO == "FL_178|FL_169|FL_133"] <- "job"
data$order_exp[data$FL_134_DO == "FL_169|FL_133|FL_178" | data$FL_134_DO == "FL_169|FL_178|FL_133"] <- "charts"
table(data$order_exp)


## order of DVs in dose response experiment

table(data$exp1_DVs_DO)
data$dose_ability_first <- I(data$exp1_DVs_DO %in% "ability_dept|rep_women")
table(data$dose_ability_first)




### ANALYSIS



## Table B9 -------------------------------------------

# Test for order effects in experiment order
data$order_exp <- factor(data$order_exp, levels = c("dose", "job", "charts"))
data$order_exp2 <- factor(data$order_exp, levels = c("job", "dose", "charts"))

summary(m1 <- lm(e1_rep_women ~ e1_gender * order_exp, data = data))
temp <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
m1$se <- temp[, "Std. Error"]
temp

summary(m2 <- lm(e1_agency_confidence ~ e1_gender * order_exp, data = data))
temp2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
m2$se <- temp2[, "Std. Error"]
temp2



stargazer(
  m1,
  m2,
  dep.var.labels = c("Represent Women", "Fulfill Mission"),
  title = "\\textbf{Pilot Experiment: Robustness Check for Order Effects.} Effects of Gender Treatment in Dose Response Experiment conditional on order in which experiment was presented in survey using Prolific sample.",
  se = list(m1$se, m2$se),
  label = "table:pilot_exp_order",
  covariate.labels = c(
    "condition: women",
    "nominee first",
    "visualization first",
    "women * nominee first",
    "women * visualization first",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)


## Table B12 -------------------------------------------
### test for effect of order of DVs
summary(m1 <- lm(e1_rep_women ~ e1_gender * dose_ability_first, data = data))
temp <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
temp
m1$se <- temp[, "Std. Error"]


summary(m2 <- lm(e1_agency_confidence ~ e1_gender * dose_ability_first, data = data))
temp2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
temp2
m2$se <- temp2[, "Std. Error"]

stargazer(m1, m2,
  column.labels = c("Represent Women", "Agency Fulfill Mission"),
  dep.var.labels = c("", ""),
  covariate.labels = c(
    "condition: women", "fulfill mission asked first",
    "women*fulfill mission asked first"
  ),
  title = "\\textbf{Pilot Experiment: Effect of Order of DVs on Treatment Effects.}",
  omit.stat = c("f", "rsq", "adj.rsq", "ser"),
  se = list(m1$se, m2$se),
  label = "table:dv_order_pilot"
)


#### Dose Response Analysis


## Figure B15 -------------------------------------------


# REP WOMEN

res2 <- dose_means_prolific(
  data = data, dv = "e1_rep_women", color_women = "#D55E00",
  color_men = "#0072B2", shape_women = 19, shape_men = 17
)
head(res2)

res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]

# a_dose_rep_women_prolific.pdf
pdf(file = "Accepted Paper Replication Archive/figures/figure_b15.pdf")

par(mar = c(6, 4, 4, 4))
plot(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], ylim = c(min(res2$lb) - .03, max(res2$ub)), axes = F,
  xlab = "", ylab = "Represent Women's Interests (0-1 Scale)",
  col = res2$color[res2$gender == "women"]
)

segments(
  x0 = 1:7, x1 = 1:7, y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1, res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1, x1 = c(1:7) + .1, y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(1, at = c(1, 3, 5, 7), labels = rev(c(
  "20%Men\n80%Women", "40%Men\n60%Women",
  "60%Men\n40%Women",
  "80%Men\n20%Women"
)), cex.axis = .8, las = 1, tck = -.02)
axis(1, at = c(2, 4, 6), tck = .02, labels = rep("", 3))

text(
  x = c(2, 4, 6), y = c(min(c(res2$lb)) - .025),
  rev(c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)

legend("topleft",
  pch = c(19, 17), lty = 1, legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ), cex = 1
)



dev.off()

### RUN AGENCY WELL

## Figure B16 -------------------------------------------


res2 <- dose_means_prolific(
  data = data, dv = "e1_agency_confidence", color_women = "#D55E00",
  color_men = "#0072B2", shape_women = 19, shape_men = 17
)
head(res2)

res2$order <- 1:nrow(res2)
res2 <- res2[rev(res2$order), ]

# a_dose_ability_prolific.pdf
pdf(file = "Accepted Paper Replication Archive/figures/figure_b16.pdf")

par(mar = c(6, 4, 4, 4))
plot(1:7, res2$coef[res2$gender == "women"],
  pch = res2$pch[res2$gender == "women"], ylim = c(min(res2$lb) - .03, max(res2$ub)), axes = F,
  xlab = "", ylab = "Agency Will Fulfill Mission (0-1 Scale)",
  col = res2$color[res2$gender == "women"]
)

segments(
  x0 = 1:7, x1 = 1:7, y0 = res2$lb[res2$gender == "women"],
  y1 = res2$ub[res2$gender == "women"],
  col = res2$color[res2$gender == "women"]
)
points(c(1:7) + .1, res2$coef[res2$gender == "men"],
  pch = res2$pch[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
segments(
  x0 = c(1:7) + .1, x1 = c(1:7) + .1, y0 = res2$lb[res2$gender == "men"],
  y1 = res2$ub[res2$gender == "men"],
  col = res2$color[res2$gender == "men"]
)
axis(1, at = c(1, 3, 5, 7), labels = rev(c(
  "20%Men\n80%Women", "40%Men\n60%Women",
  "60%Men\n40%Women",
  "80%Men\n20%Women"
)), cex.axis = .8, las = 1, tck = -.02)
axis(1, at = c(2, 4, 6), tck = .02, labels = rep("", 3))

text(
  x = c(2, 4, 6), y = c(min(c(res2$lb)) - .025),
  rev(c("30%Men\n70%Women", "50%Men\n50%Women", "70%Men\n30%Women")),
  cex = .8
)

axis(2, at = seq(0, 1, .05), las = 2)

legend("topleft",
  pch = c(19, 17), lty = 1, legend = c(
    "% women treatment",
    "% men treatment"
  ),
  col = c(
    res2$color[res2$gender == "women"][1],
    res2$color[res2$gender == "men"][1]
  ), cex = 1
)



dev.off()


#### NOMIINEE ANALYSIS

## Table B17, pilot results -------------------------------------------


class(data$party2)
data$party2 <- factor(data$party2, levels = c("rep", "ind", "dem"))
levels(data$party2)

class(data$e2_gender)
data$e2_gender <- factor(data$e2_gender, levels = c("male", "female"))
levels(data$e2_gender)

m1 <- lm(e2_nominee_confidence ~ e2_gender, data = data)
temp <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
temp
m1$se <- temp[, "Std. Error"]

m2 <- lm(e2_nominee_confidence ~ e2_gender * party2, data = data)
temp <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
temp
m2$se <- temp[, "Std. Error"]

m3 <- lm(e2_nominee_confidence ~ e2_gender * female, data = data)
temp <- coeftest(m3, vcov = vcovHC(m3, type = "HC1"))
temp
m3$se <- temp[, "Std. Error"]

m4 <- lm(e2_pres_confidence ~ e2_gender, data = data)
temp <- coeftest(m4, vcov = vcovHC(m4, type = "HC1"))
temp
m4$se <- temp[, "Std. Error"]

m5 <- lm(e2_pres_confidence ~ e2_gender * party2, data = data)
temp <- coeftest(m5, vcov = vcovHC(m5, type = "HC1"))
temp
m5$se <- temp[, "Std. Error"]

m6 <- lm(e2_pres_confidence ~ e2_gender * female, data = data)
temp <- coeftest(m6, vcov = vcovHC(m6, type = "HC1"))
temp
m6$se <- temp[, "Std. Error"]



stargazer(m1, m2, m3, m4, m5, m6,
  dep.var.labels = c(
    "Nom\\ Confidence",
    "Pres.\\ Confidence"
  ),
  title = "\\textbf{Pilot Study: Effects of Gender Treatment in Nominee Experiment.} The table below shows
          average treatment effects of using female pronouns relative to male in the nominee experiment pooled
          across all doses. Models (a) and (data) show average effects in the entire sample.
          Models (b)-(c) and (e)-(f) condition on respondent party and gender, respectively.",
  se = list(m1$se, m2$se, m3$se, m4$se, m5$se, m6$se),
  label = "table:nominee_pilot",
  covariate.labels = c(
    "condition: female pronouns",
    "Independent", "Democrat", "female pronouns x Independent",
    "female pronouns * Democrat", "female respondent",
    "female pronouns * female respondent",
    "(Intercept)"
  ),
  omit.stat = c("f", "rsq", "adj.rsq", "ser")
)

